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Abstract 

Bilayer lipid membranes (BLMs) are an essential component of all biological systems, forming a 
functional barrier for cells and organelles from the surrounding environment. The lipid molecules 
that form membranes contain both permanent and induced dipoles, and an electric field can induce 
the formation of pores when the transverse field is sufficiently strong (electroporation). Here, a 
phenomenological free energy is constructed to model the response of a BLM to a transverse static 
electric field. The model contains a continuum description of the membrane dipoles and a coupling 
between the headgroup dipoles and the membrane tilt. The membrane is found to become unstable 



through buckling modes, which are weakly coupled to thickness fluctuations in the membrane. The 
thickness fluctuations, along with the increase in interfacial area produced by membrane buckling, 



in 

increase the probability of localised membrane breakdown, which may lead to pore formation. 



The instability is found to depend strongly on the strength of the coupling between the dipolar 
headgroups and the membrane tilt as well as the degree of dipolar ordering in the membrane. 
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I. INTRODUCTION 

When amphiphilic lipid molecules are dissolved in solution, the molecules can self- 
assemble into a bilayer structure with the hydrophilic headgroups of the molecule shielding 
the hydrophobic hydrocarbon tails from the surrounding water. Many biological processes 
and components that occur in the cell depend on the membrane: membrane-bound proteins, 
endo/exocytosis, lipid rafts and ion channels are just a few examples (l|. The many different 
species of lipid found in a cell membrane share the same general structure: a polar headgroup 
attached to a non-polar hydrocarbon tail region. When a bilayer lipid membrane (BLM) 
forms, the non-polar tails make up the core of the membrane, with the dipolar headgroups 
forming the membrane surface [2| . Both parts of the molecule react to electric fields. When a 
strong electric field is applied transversely across a membrane, reversible electric breakdown 
can occur. The breakdown is characterised by an increase in the measured conductivity 
due to the rapid increase in the transit of ions across the membrane 3L |4j . This increase 
in permeability is attributed to the development of transbilayer pores 5] , which may close 
upon removal of the electric field, allowing the membrane to recover. This phenomenon has 

nil 

been termed electroporation |6|, |7|. 

The theoretical work on electroporation and electrical breakdown can be viewed as be- 
longing to two distinct branches: one approach uses the Smoluchowski equation to describe 
the evolution of a distribution of pores with an assumed energy for pore formation |8l-ll6l|. A 
density of pores are modelled drifting through radius space as function of time, generated by 
a source term including the effect of the field. This method has been successful at predicting 
pore radii, lifetimes and densities, but does not model the mechanism of pore formation [7|. 
A different approach is required to understand how pores form and what membrane proper- 
ties inform this process. The simplest approach is to coarse-grain the BLM to a continuum 
membrane driven unstable by an electric field. Early work by Crowley [17] modelled the 
hydrocarbon core of the membrane as a dielectric slab with finite shear modulus and finite 
elastic compressibility, but estimates a critical transmembrane voltage an order of magnitude 



larger than experimental values |7|. The model of Lewis 



a. 



also models the membrane as 



a dielectric slab, but includes a Maxwell stress tensor which relates the dielectric constant 
to strain in the membrane, however also finds a critical transmembrane voltage larger than 
those experimentally reported [6|, similar to Crowley. 



These models neglect the fluid bilayer structure of the membrane, and thus neglect im- 
portant mechanical properties such as a vanishing shear modulus and the bending rigidity. 
The Helfrich-Canham Hamiltonian and its variants [2| are frequently used to model con- 



formational changes to the membrane. Sens and Isambert [19[ adapted these methods and 
considered the minimisation of the difference between the stressed and unstressed areas of a 
membrane in an electric field. The authors imposed a force from the electric field on an un- 
dulating membrane and calculated the unstable undulatory wavelength and corresponding 
growth rate, although the model used neglects any thickness variation in the membrane. The 
stochastic thermal undulations proposed as a mechanism for pore formation by Movileanu 
et al. 20( are hindered by a large energy barrier (91/cbT) and neglect the effect of the field 
on the membrane. Membranes are only nm in thickness and the process of membrane break- 
down under electric fields occurs over short time scales, which makes experimental study of 



pore formation difficult. 

Molecular dynamics (MD) studies o 



electrical behaviour of membranes 
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' membranes have been used extensively to study the 



-24|. These simulations provide molecular- level detail 
on a picosecond time scale. However MD can only simulate a very small area of membrane 
for a short time. The transbilayer pores opened during electroporation can last for up to 



[G 



ms [6[ before closing, which MD simulations cannot capture. Experimental studies range 
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261 ] to conductivity measurements 



from measurement of the transmembrane current 

using salt-filled vesicles [27[. Recent developments in video microscopy and fluorescence 
lave enabled the direct visualisation of giant unilamellar vesicles exposed to an electric field 
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-|30|. in which pores can be directly observed. 

In this work, we develop a comprehensive, mesoscopic analytical approach, which in- 
cludes a mechanical coupling between the orientation of the dipole on the surface and the 
membrane surface tilt. This should destabilise the membrane as the headgroups seek to 
align with the field, rather than shield the hydrophobic core of the membrane from the sur- 
rounding fluid, which is their equilibrium position. As the headgroups tilt, the membrane 
will tilt to try to restore the equilibrium position, which can introduce an instability in the 
membrane not noted previously in the literature. We study this instability by performing a 
linear stability analysis of the free energy. This perturbative approach will not capture the 
inherently discontinuous process of pore formation, but will predict the onset of instability 
in the membrane. The instability occurs through deformational modes involving thickness 



fluctuations in the membrane, which increases the probability of localised breakdown and 
therefore of pore formation. Applying a field to the membrane breaks the up/down symme- 
try of the membrane therefore it is important we include a description of the bilayer which 
allows each monolayer to be independently deformable. 

In Section [TT] we construct the free energy including terms associated with mechanical 
deformation and introduce the description of the dipolar headgroups. Section IIIII presents 
the qualitative analysis of the model, Section llVl presents results from numerical calculations, 
and we conclude the paper and discuss possible future work in Section [V] 

II. MODEL 

A. Geometry 

We consider a planar bilayer lipid membrane suspended horizontally in water with an 
electric field applied such that the field is perpendicular to the unperturbed membrane 
surface. The membrane is modelled as a dielectric, fluid membrane at zero tension with 
a non-zero area stretching modulus. The generalised three-dimensional (3D) free energy is 
given in Appendix |Aj but to illustrate the basic principal and obtain analytically tractable 
solutions we assume a one- dimensional modulation in the x direction (Figured]). 
Here h± denote the positions of the upper (+) and lower (-) membrane surfaces, t± is 
the thickness of the upper and lower membrane leaflets, to is the unperturbed monolayer 
thickness, and s is the displacement of the dividing surface between the monolayers. 

B. Conventional free energy 

We construct a phenomenological free energy per unit area. The first contribution is the 
energy associated with mechanical deformation of the membrane; 
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Here ka is the area compressibility modulus, Kb is the bending rigidity, 7 is the surface 
tension and t is the initial leaflet half-thickness. The primes represent differentiation with 



respect to the x direction. These terms are equivalent to those used by Huang in [3jJ, and 
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FIG. 1: (a) An unperturbed bilayer, (b) the bilayer following a deformation. 



have been adapted from the Helfrich-Canham Hamiltonian. The surface tension in our model 
is not equivalent to a frame tension, which acts in the bilayer midplane. Instead, it restricts 
variation in interfacial area of each leaflet separately and hence can describe peristaltic 
deformations. The area compressibility term allows the two monolayers that form the bilayer 
to be independently deformable [32] . This has a strong effect on the relaxational dynamics 
of the membrane, but in the static case considered here these deformations will equilibrate 
effectively instantaneously, meaning we can minimise over s at this stage without loss of 
generality. [Note that the large bending modulus (i% b rs 10 k^T) precludes renormalisation 
of these elastic constants [33|.] 

We also include the dielectric energy f d : 



Qi 



f d = -± [e m E 2 m (h+ - h.) + e w E\ (L - h+) + e w El (h_ - L 



(2) 



Here eo is the dielectric constant of the vacuum, E m is the field in the membrane, and e m 
and e w are the dielectric constants of the membrane and water respectively. E± is the field 
at the upper (+) and lower (-) membrane surface. L is the upper and lower limit of the 
system. The form of Eq. ([2]) implies the field will cause a uniform compressive force on 
the membrane. This 'electrostrictive' force has been found to have a quadratic dependence 
on the applied transmembrane voltage and has a small effect (~ 1% fractional thickness 



change) for the voltages used here [34|-|36|. 



C. Dipolar Headgroups 

The dipolar headgroups of the lipids are defined by a three dimensional vector p; 

p = p z + m 

= pcos(8) e z + m, (3) 

where p = |p|, 9 is the angle the dipole makes with the z direction, e z is the unit vector in 
the z direction and m is a vector representing the in-plane dipole moment. The value used 
for p is the effective magnitude of the headgroup dipole moment, calculated by Raudino and 
Mauzerall 37j, which includes the screened charges and conformation of the headgroup. 

In an unperturbed bilayer lipid membrane, the headgroup of each phospholipid molecule 
lies at an average angle of 9$ to the membrane normal, hinged about the uppermost carbon 



atom [38( . This natural tilt of the headgroup arises from a balance between the dipole-dipole 
interactions, the shape of the molecule and the need to shield the non-polar hydrocarbon 
chains from the water. We perturb about the equilibrium position: 

6± = 6 ± + 86± 

where 

#0+ = #o 

O - = 7T ± (A and B) . 

The dipolar orientations between the two leaflets are weakly coupled by the Coulomb 
interaction, which prefers antiparallel orientations. In principle this degree of freedom allows 
for rich phase behaviour and dynamics. Since the dipole orientation is coupled to the applied 
field, the relative orientation of the dipoles on the upper and lower membrane leaflets will 
affect the membrane behaviour under an electric field. Here we consider only parallel and 
antiparallel orientations. We refer to the case of the dipoles pointing in opposite directions as 
antisymmetric (A), and the case where the dipoles point in the same direction as symmetric 
(B), as shown in Figure |2j 




FIG. 2: Symmetric Sz antisymmetric dipolar orientation between leaflets. 

Tilting the dipole relative to the membrane surface will cost energy, reflected in the 
following free energy per dipole; 
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2 
Here /t p is the dipole-membrane coupling modulus and e x is the unit vector in the x direction. 
This term is inspired in part by the work of Lubensky, Chen and MacKintosh in 
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411, 



and explicitly penalizes change in dipole orientation relative to the membrane surface tilt. 

So far p and m have referred to a single dipole. To treat larger membrane areas, we 
extend these to continuum variables: p becomes a dipole moment density p (dipole moment 
per unit area, with dimensions Cm -1 ). We coarse-grain m into m = (m), the average 
orientation of dipoles within a small area. To distinguish between changes in orientation of 
the dipoles and changes in alignment within the coarse-grained area, we separate m into 



two components; a unit vector rh = m/ | m| , representing the average orientation and an 
amplitude rh = |m|, which represents the degree with which dipoles within the given area 
point along rh. If all the dipoles within the area point along rh, then rh — 1. If rh — 
then the dipoles within the area are completely disordered. As we are considering a one- 
dimensional modulation in the x direction we can set rh = e x and allow rh to vary. We write 
rh as an equilibrium value rh , and the deviation 5rh±, from this equilibrium value. 

rh± = rh + 5rh±. 

The equilibrium value m arises from a competition between the dipole-dipole interactions 
and their thermal fluctuations. We allow 5fh± to vary independently between the leaflets. 

Changes in dipole alignment are penalised by a susceptibility Xm, leading to the free 
energy density 

fx = ^Sfhl (5) 



A three-dimensional construction of f x is given in Appendix |A] As in Andelman et at. 42] 
we include a term coupling the dipole alignment with the surface curvature: 



f c = -^[{h'i)5m' + +(h'L)5m'_] 



(6) 



where 7 C is the relevant modulus. 

The final contribution to the free energy, fg is a free energy of two uncoupled dipoles in 
an electric field. The free energy per unit area is 
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antisymmetric case (A) 



(7) 



601) + sin (6» ) (58+ + 56L) symmetric case (B) 



where we have assumed equal fields in the upper and lower water regions, E+ = EL = E. 
Here we have also assumed that the field acting on the dipoles is equivalent to the field 
at the membrane surface, not the field in the membrane core, which differs by two orders 
of magnitude. The dipoles are in an aqueous environment which is very different from the 
hydrocarbon core of the membrane 18(. The relative dielectric constant of bulk water is 
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?s 80, whereas for the membrane interior it is ~ 3. The field is assumed to only act between 
the two membrane surfaces hence the effect of ionic screening can be neglected. Here 9± 
refers to an average tilt angle within the region of coarse graining. 

Now that the free energy has been constructed, it can be subjected to a linear stability 
analysis. This will predict the onset of a static instability in the membrane. In an ex- 
perimental system, the instability is complicated by the dynamics of the surrounding fluid 
(likely to contain many ions, particularly near the membrane surface) and by the dynamics 
of the bilayer itself, which behaves as a two-dimensional fluid. To capture the dynamical 
behaviour of the instability predicted by our model, we would need to include both the 



hydrodynamic flows of the fluid and membrane 32|, |43| and the movement of charges in the 



solution 
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451 ] . These are both non-trivial extensions due to the explicit definition of the 



dipoles in our model Eq. (jM7|) and hence beyond the scope of this paper. 

III. QUALITATIVE ANALYSIS 

A. Fluctuation free energy 

After constructing the free energy we change variables to modes that characterize the 
bilayer as a whole: 

u = h+ — h- peristaltic mode 

h = h + + /i_ bilayer mode 

A = 59 + — 59^ difference in dipole tilts 

X = 59 + + 59 '_ mean dipole tilt. 



From Equations (tTJlTj), the free energy can then be expressed as; 

/ = — \u + h J + - \u + h J + — (Sm + + Sm_) 



K A (u 2 2 \ e e w £ 2 /e w 
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^ [(/*" + u") 6m' + + (fc" - u") 6m'_] 



(8) 
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+p E I cos 2 °^ A S + sin (0 O ) S J antisymmetric 

+p E ( cos 2 si A S + sin (0 O ) A J symmetric. 

The system has six remaining degrees of freedom; 6m±, A, E, w and h. We minimise / 
over the dipolar tilts A & E. The minimised values of these tilts are given in Appendix 
To quadratic order, the resulting free energy for the symmetric and antisymmetric case is 
identical. Figure [3] shows the dipole configurations for the bilayer & peristaltic modes. The 
variables u and h are expanded as small perturbations about the flat state; 

u — Uq + Su h = h + Sh (9) 

where 

d x u = d x ho = 0. (10) 

The equilibrium value of the peristaltic mode uq includes the electrostrictive thinning implied 
by Eq. ©. 

The free energy will then consist of fo (uo,ho) and the perturbation 5f (5u,Sh). The 
flat state fo (uo, ho) describes the membrane after the application of an electric field in the 
absence of undulations. We expand the perturbation to second order in 8h, 5u and Sm±. 
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FIG. 3: The dipole orientations in an applied field. The field applied to each membrane ((f) = 0.03) 
is equivalent to a voltage drop of 0.07V across the membrane. The parameters used are given in 
Table [III The dipole orientations were generated using the equations for the minimums of A and E 
found in Appendix O The functional form of the bilayer deformation modes u and h are imposed 
in order to display pure modes. 
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The stability of the system is determined by 5f; 

Sf = * (^h" 2 + 5u" 2 ) + 1 (5u> 2 + 5h' 2 ) + ^ (^ + 5m 2 _) 



^/^\ + *^f Mm 
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[(5/i" + 5m") 5m! + + (5/i" — 5m") Sm'_] . 



when 5f < the system becomes unstable. 



B. Fourier Expansion 



Upon Su and 5/i decomposing into Fourier modes, 



5u = J2 u q t e-^ x/t0 



8h = J2 h q t e-^ x/t0 



22 s qtoe 
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(12) 
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where q is a normalised Fourier wavenumber given by q = qt , The free energy fluctuation 
is given by 



where 
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[(h* q + u* q )6m+ + (h* q -u* q )Sm-]. 
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and 

_ k p m 1 _ 7t 2 _ _ ^a^q 

^p — o s — — — a a — 

K b ^ K b 

Xratl pE COS (gg) tg 

The dimensionless parameters reference all energies to the membrane bending energy K b . 
The dimensionless potential <fi compares the strength of the field on the dipoles with the 
membrane bending rigidity. For </> <C 1, we would therefore expect the system to be stable 
as the bending rigidity will be much stronger than the field on the dipoles. We predict that 
the system will only become unstable for <ft w 0(1) where the strength of the field on the 
dipoles can overcome the bending rigidity. The effect of the field through <fi is regulated by 
a p ; the two terms appear concurrently in Eq. (I14p . This is because a p controls the way the 
membrane 'feels' the field through the dipole- membrane coupling Eq. (J3J). 



C. Matrix Representation 

We rearrange Eq. ( Tl4|) in the quadratic form 

S fq = o ^ ' M< ? ' V Where V = ( U 9' ^«> 5 ™t' S ™q) ^ 

where v is complex. The form and components of M g are given in Appendix [B] We calculate 
the eigenvalues X q and eigenvectors e\ of M ? , which describe the eigenmodes of the system. 
The membrane is unstable when any of the eigenvalues become negative. The associated 
eigenvectors e\ determine how much the bilayer (h q ), peristaltic (u q ) and dipole alignment 
(Srri^) modes are involved in each eigenmode. If the eigenvector is dominated by one of 
these pure modes, we will refer to the eigenvalue as distinctly associated with the pure mode 
that dominates its eigenmode. The matrix M g can also be used to calculate the fluctuation 
spectrum of the modes; 

{v i v j ) = a(M- l )..k B T, (16) 

where a = 2 tt 2 /L 2 is a constant related to the Fourier expansion. 

D. Analytical Features 

We can gain some qualitative information about the membrane stability by analysing 
the fluctuation free energy, Eq. ( TT3|) . The dimensionless parameters a a, cr s , a p , a x and a c 
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represent the effects of the area compressibility, surface tension, dipole-membrane coupling, 
dipole alignment and dipole alignment-membrane coupling on the free energy, scaled by the 
bending energy. The typical values for these ratios, calculated using the parameters given in 
Table HI are given in Table HT1 Since a a > 1, the membrane is more susceptible to bending 
than compression. The ratio o p controls how the membrane 'feels' the applied electric field 
through the dipoles, as it contains both the dipole-membrane coupling modulus k p and the 
average dipole alignment rho- Since this ratio a p is typically < 1, this effect is dominated 
by membrane bending. The ratio containing the surface tension, a s , which using the values 
of the parameters in Table [I] is < 1. This will cause the membrane to bend with surface 
gradients rather than surface curvature. 

The model should be stable (5f q > 0) for a flat membrane (q — 0). For zero field {(f) = 0) 
and in the limit q — > 0, the thickness variations (peristaltic modes u q ) are penalised by the 
area compressibility a a, while the bilayer modes vanish. Since u q is stabilised against small 
changes in q, we expect the instability to progress via bilayer (hgj modes. If we consider 
only the bilayer modes, the instability occurs when the term proportional to q 2 , which can 
be considered to be an effective surface tension, changes sign: 

(17) 



1 + a p /a s 



The critical potential C is therefore approximately proportional to the ratio a p [53| . The 
ratio a p contains two important independent parameters: /%,, the strength of the dipole- 
membrane coupling and the average dipole alignment m . The membrane is stabilised upon 
increasing either k p or rho. Increasing k p 'stiffens' the dipoles against movement away from 
their equilibrium position, while increasing rho increases the proportion of dipoles within a 
membrane patch that are aligned along the x axis and therefore constrained by Eq. ([6]). As 
reducing m lowers C , we can infer that an instability is more likely to occur in a membrane 
region in which the dipoles are disordered (i.e. rh is smaller). The generalised 3D version 
of this term considers the orientation of the dipole in the complete plane, rather than just 
the x-direction, which could describe two dimensional modulations. We expect this to be 



energetically more costly at the onset of instability [39]. 

For the values used given in Table HT1 (a* = 0.20, a p = 0.23), 4> c ~ 0.15 which is equivalent 
to a voltage drop of 0.22V across the membrane. This simple qualitative estimate gives 
values for the critical potential that are close to the experimentally reported range of (0.2- 
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IV) |6j. The inclusion of the coupling between the bilayer and peristaltic modes should 
raise the critical potential by transferring energy from the bilayer modes into the peristaltic 
modes. 

IV. NUMERICAL CALCULATIONS 

A. Parameters 

All the parameters used to obtain these results are given in Table [H The upper portion 
of the table contains the parameters that have been obtained from experimental studies, 
whereas those in the lower portion are not accessible by current experiments. Those pa- 
rameters that cannot be experimentally measured have either been estimated from physical 
reasoning (% m , 7 C and m ) or extrapolated from simulation results (k p ). Of these constants, 
we expect only k p to have a significant effect on the stability of the membrane due it being 



present in Eq. [TO T 



re estimate of k p comes from the distribution of dipole angles obtained 



by Bockmann et al. 38[. The authors perform MD simulations of lipid bilayers and measure 
the angles formed between the dipole of the lipid headgroup and the surface normal. This 
result has been produced in other studies, using different simulation methodologies |46|, |47|. 
We assume the width of this distribution is governed by one degree of freedom and then cal- 
culate the stiffness with which the dipole hinges around the equilibrium position, assuming 
it bends as a Hookean spring. 

B. Eigenvalue stability 

Figure H] shows the eigenvalue associated with the bilayer modes h q . As the rescaled 
potential <fi is increased, the eigenvalue A* decreases. For <fi = 0.21 the eigenvalue has a 
significant negative portion indicating that has passed through the point at which the 
membrane first becomes unstable. For this unstable value of <p the instability begins at 
q = and reaches a lowest value at q = 0.6, which corresponds to peak to peak spacing 
of 26nm. The eigenvalue associated with the peristaltic undulations u q behaves identically 
to the eigenvalue for h q , but has a ^-intercept of 2a a, hence this branch will never become 
unstable for q G [0,1]. The membrane can therefore become unstable only through the 
bilayer modes of undulation, as suggested in the previous section. While this is not an 
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TABLE I: Parameters used for calculations 



Symbol 




Name 


Value 


Ref 


Lipid 


KA 




Area compressibility 


0.f4Jm" 2 


flj 


POPC 


«6 




Bending rigidity 


0.4 x 10~ 19 J 


[48] 


DMPC 


7 




Surface tension 


1.5 x lCr 3 Jm- 2 


[31] 


GMO 


0o 




Equilibrium dipole orientation 


60° 


[38] 


POPC 


P 




Dipole moment per unit area 


1.1 x lO^Cm -1 


[18] 


POPC 


to 




Monolayer thickness 


2.5 x l(T 9 m 


py 


POPC 


ftp 




Dipole-membrane coupling" 


1.5 x lCr 2 Jm- 2 


[38] 


POPC 


Xrn 




Dipole alignment susceptibility 


3.45 x l(r 2 Jm- 2 


- 


- 


7c 




Dipole alignment-membrane coupling 


0.0-0.6 x 10~ 19 J 


- 


- 


fh 




Average degree of dipole alignment^ 


0.3-0.4 


- 


- 


"The di 


po 


e-membranc coupling modulus, k p , was generated from the distribution of hcadgroup angles 


given in 


38 


. The width of this distribution was assumed to be to be governed by one degree 


of freedom. 



This estimate is then multiplied by the number of lipids per unit area, nu v (~ 10 19 lipids/m 2 ) [l[. 
6 Thc dipole alignment susceptibility, Xm, was estimated as \ m ~ kbTriup, due to the entropic cost of 

constraining the dipole alignment 
The functional form of a c compares the dipole alignment-membrane coupling strength 7 C with the bending 

rigidity Kb- As the effect that Eq. ([6]) regulates is not seen experimentally, we can assume that "f c < Kb- 

However the range of "f c is chosen in order to provide examples where "f c > Kb as well as the physically 

expected range. 
d The range of the dipolar alignment ffiQ used was judged to be a reasonable estimate of the degree of dipole 

alignment in a membrane. 



obvious route to transmembrane pore formation, bilayer modes have been observed to play 

|, as well as occurring in the 



an important role in MD simulations of electroporation 
theoretical model of Sens and Isambert 19]. 

The locus of instability as a function of q and <fi is shown for two values of a p in Figure[5]^a). 



For ij p = 0.23 the membrane becomes unstable at the critical potential of 



0.16. This 



corresponds to a voltage drop across the membrane of roughly 0.24V, which is in the range of 
values (0.2-1V) for the onset of electroporation seen in both experiments 
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and simulations 
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FIG. 4: The eigenvalue A^ associated with the bilayer undulations h q as a function of q. The 
shaded region is unstable. 



221 ] . This is slightly higher than the qualitative estimate obtained in the previous section. 
The instability in our model is likely to be relieved by a change of state of the membrane. The 
formation of transmembrane pores can achieve this by allowing ions to permeate through 
the system, reducing the electric field across the membrane. 

Increasing the dipole-membrane coupling strength a p to 0.34 increases the critical poten- 
tial to <p c = 0.22 (0.34V). This is again slightly larger than the value predicted by Eq. (EC 
(<f) c = 0.20), but is also smaller than predicted by a linear increase of <f) c with a p . 
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C. Dipole alignment- membrane coupling 

The effect of the dipole alignment-membrane coupling a c on the membrane stability is 
shown in Figure [5(b). Increasing a c does not affect the onset of the instability at long 
undulatory wavelengths (q = 0), where the cubic dependence on the wavelength (g 3 ) of a c 
in the free energy is outweighed by the quadratic dependence (q 2 ) of the terms which cause 
the instability. The variation of the dipole alignment-membrane coupling a c does affect the 
shorter wavelength (q — >■ 1) structure of the instability, where the cubic and quadratic terms 
become comparable. This will affect the behaviour of the membrane if a field larger than the 
critical value is applied rapidly. Overall 7 C only weakly affects the stability of the membrane, 
over the range 7 C = — 1.5 ac&. 

Figure[5(b) shows the effect of the dipole alignment-membrane coupling a c on the stability 
of the membrane and thus only the effect of a c on the bilayer modes. The variation of a c 
affects the peristaltic (u q ) modes differently, as shown in Figure [61 The peristaltic modes are 
stabilised at higher q by an increase in the dipole alignment-membrane coupling a c , whereas 
the bilayer modes are destabilised. 

The dipole alignment-membrane coupling energy (Eq. [6]) couples the membrane deforma- 
tion modes (h q and u q ) with the dipole alignment modes Smi^. The eigenvalues A^, whose 
eigenvectors are dominated by the dipole alignment modes Srri^ are shown in the lower sec- 
tion of Figure [61 The q = limit of these eigenvalues is governed by and proportional to 
<t x , and the eigenvalues only deviate from this value at larger q. Significant change in the 
eigenvalues can only be seen for values of o c > 1 . These large values of a c are unlikely to be 
physically realisable as they require 7 C ~ O (ft fc ). The magnitude of 7 C cannot be measured 
directly by current experiment. However as the degree of dipole alignment is observed to 
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49| it is fair to assume 7 C <C k&. 



have little correlation with membrane bending 

For the 5m + - and <5m_-dominated eigenvectors, the corresponding eigenvalues A^ are 
degenerate at q = and then separate as q increases. Increasing the strength of the dipole 
alignment-membrane coupling a c increases the amount by which these eigenvalues deviate. 
The corresponding eigenvectors are not associated with Sm + and Sm_ independently, but 
rather with both modes equally; however the more stable eigenvalue has a slight contribution 
from the bilayer mode h q , whereas the less stable eigenvalue has a slight contribution from 
the peristaltic modes u q . The difference between these modes explains why the eigenvalues 
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FIG. 5: The locus of instability as a function of q and <j> for various values of (a) the dipole- 
membrane coupling a p and (b) the dipole alignment-membrane coupling a c . The shaded region is 
unstable. 

associated with the bilayer and peristaltic modes react differently to increases in a c 

D. Fluctuation spectrum 

Whereas the variation of a c has a small effect on the overall stability of the membrane, the 
interaction of a v with the fluctuation modes provides a test of the model. The fluctuation 
spectrum of the modes can be calculated using Eq. (TT6l) . As expected, the fluctuations of 
the bilayer modes are much larger in magnitude than the peristaltic modes. This is because 
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FIG. 6: The eigenvalues for both the peristaltic and bilayer modes (upper panel) and the dipole 
modes (lower panel) as functions of q. As a c is increased at constant (j) (=0.21), the bilayer and 
peristaltic modes react differently. The lower branch of the bilayer modes becomes stable for q > 1. 
The dipolar modes show little difference for field strengths above or below the critical field strength 
{(p c = 0.16), but respond strongly to changes in the dipole alignment-membrane coupling a p . As q 
is increased, the dipolar modes deviate from their initial value, the direction determined by their 
weak association with either the bilayer modes h q or the peristaltic modes u q . 

the peristaltic modes are dominated by the strong stretching modulus, a a- For small a a the 
fluctuations of the peristaltic modes grow to match the fluctuations of the bilayer modes. 
For zero applied field {\h q \ 2 ') ~ 1/ (g 4 + a s q 2 ), as expected for a flat membrane [2[. The 
ratio of fluctuations of the bilayer modes at = 0.15 to = 0.0, (\h q \ 2 ) , lg / (\h q \ 2 ) , =00 is 
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TABLE II: Dimensionless ratios according to measured data from Table U 



Symbol 


Formula 


Value 


o-.l 


Kb 


21.3 


08 


i4 

K b 


0.2 


°V 


Kp ?71q Lq 
Kb 


0.23-0.35 


CT \ 


Kb 


5.4 


£T C 


Kb 


0-1.5 





pE cos(0o)*o 

Kb 


0-0.35 



displayed in Figured for various values of the dipole-membrane coupling strength a p , show- 
ing that the application of a field increases the magnitude of the fluctuations. Conversely, 
the dipole-membrane coupling a p 'stiffens' the membrane, and reduces the fluctuations as 
reflected in Figured 

E. Eigenvector composition 

Figure [8] shows the unstable eigenvalue and the associated normalised eigenvector. For 
zero applied field, e\-e h = 1 and e\-e u = 0, where e h and e u are the unit vectors representing 
the pure modes h q and u q respectively, hence the eigenvalue is associated only with the 
bilayer modes. For <ft > <fr c the amplitudes of the eigenvector components vary with q 
and the contribution from the peristaltic mode e A • e" increases due to the coupling term 
present in Eq. ( 1T4|) . This term has a similar field dependence to the effective surface tension 
term present in Eq. ( jl~4"|) which induces the instability, so the eigenvector mixing increases 
dramatically for > <f) c . The change in the eigenvector components is small compared with 
the initial composition, therefore we can consider the eigenvalue as distinctly associated the 
bilayer modes. This behaviour is mirrored in the eigenvalues of the pure peristaltic modes, 
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FIG. 7: The ratio between the fluctuations of the bilayer (hq) modes at <p = 0.15 to <fi = as a 
function of g. 

with the bilayer mode contribution (£\ ■ e h ) increasing slightly for <fi > (p c and increasing q. 
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FIG. 8: The unstable eigenvalue and corresponding eigenvector e\ as a function of q. e\ ■ e and 
e\ ■ e u are the contributions to £\ of the bilayer and peristaltic modes respectively. 

V. DISCUSSION & SUMMARY 

We have constructed a model of a planar membrane in an electric field, which contains an 
explicit coupling between the orientation of the dipolar lipid headgroups and the membrane 
shape, thus coupling the application of the field to the membrane shape in a way not seen 
previously in the literature. The phenomenological model contains only harmonic terms 
which are subjected to a linear stability analysis. This model becomes unstable as the 
applied field is increased, with a critical potential that matches those seen in experiment 
and simulation J7|. A simple formula (Eq. ( ITT)) ) has been found that gives a reasonable 
estimate of the critical potential related to a minimal number of model parameters, which is 
useful as decreasing the number of parameters used decreases possible sources of error. The 
instability depends strongly on m , the average alignment in a membrane patch, with the 
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instability occurring for smaller fields for disordered membranes of smaller mo- As dipole 
alignment will vary dynamically in a physical system, the membrane is more likely to become 
unstable in disordered patches. This means the model captures some of the stochastic nature 
of membrane breakdown and pore formation. The instability also depends strongly on k p , 
the strength of the coupling between the dipolar headgroups and the membrane core. This 
is likely dependent on the combination of lipids in the bilayer. Since variations in m® or k v 
have a significant effect on the critical potential <f) c , these would be good parameters with 
which to test the model. 

Because the process of membrane breakdown requires a rupture to form in the mem- 
brane, it cannot be fully modelled by any continuum theory. Despite this, the instability 
studied in this work can be linked with the formation of defects within the membrane and 
therefore the formation of transmembrane pores. From Figure [U the unstable eigenvector 
shows that the instability is dominated by the bilayer modes but approximately 2.5% of the 
instability involves the peristaltic modes. This induces a periodic thinning which destabilises 
the membrane. Evans, Waugh and Melnick J50( found using micropipette aspiration that a 
membrane can only support thickness changes of ~ 4% before rupture. To induce a frac- 
tional thickness change of this magnitude using the peristaltic undulations produced by the 
instability requires the bilayer modes to have an amplitude of 6nm. This is above the size 
that would be produced spontaneously by thermal fluctuations, but after the application of 
an electric field, the bilayer modes become unstable and this amplitude could be achieved. 
A membrane defect is then more likely to form at the troughs of the peristaltic undulations, 
where the membrane is thinnest. This defect could then go on to form a pore, or rupture the 
entire bilayer. The most unstable undulation wavelength is 26nm (Figure H]), comparable 



with the average pore-pore separation reported in 5jJ, consistent with this hypothesis. 
The parameters k p and 7 C , both unique to this model, provide opportunities to make 

predictions and test this model. An obvious extension to the model would be to allow for 
;he full rotation of the dipole distribution, which could lead to non-trivial pattern formation 
3914411 • Our calculations only calculate the static instability. To fully model the dynamical 

behaviour of the instability predicted in this work, we would need to include both the 

hydrodynamic flows of the fluid and membrane 
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solution 
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431 ] and the movement of charges in the 



45] . Coupling hydrodynamic flows to the movement of the membrane would 



be expected to push the instability to smaller wavelengths (larger q) [521 ] . 
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Appendix A: 3D free energy 

The general 3D form of the free energy of membrane deformation; 



ft 6 



U = f ((v 2 /* + r + (v 2 M J + ~ ((VM 2 + (VM 2 ) 



i 



ms- 1 -*** 



1 + £ n V 2 s 



(Al) 



For the free energy associated with the dipole surface coupling f p 



K, 



fp = y i ((p+ - p+) • ") p+ • v/i + - (p+ • ( v/l + ~ v/l o+)) Ip±+I 



+ 



((p.- $-)•*) 



p • V/i- 



: P _ • (v/i_ - v/i -)) ip_ 



where 



p ± = p - (p • n) n 



(A2) 



and the hatted variables are normalised. ho± is the initial surface gradient and p* is the 
perturbed dipole vector. 

f c is the free energy of the dipole alignment-membrane coupling ; 



7c 



fc = f ( (v 2 /^+ V • p+) + (V 2 /i_ V • p_) ) 



(A3) 



/ x is the energy punishing dipole alignment; 



fx 



Xr, 



P±+ - P±+) + (P_|__ - P±-J 



(A4) 



fd has the same functional form but is integrated over three directions instead of two. 
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Appendix B: Matrix representation 



The matrix constructed for Eq. (IT51) is given by; 



M„ 



Mu 





M13 








M 16 





-M 16 





M n 





M 13 


-M 16 





M 16 





M 13 





M 33 








M 16 





M 16 





Mia 





M 33 


-M 16 





-M 16 








-Mi 6 





-M 16 


M 55 











M 16 





Afie 








M 55 











M 16 





-Mi 6 








M 55 





-M 16 





M 16 














M 55 



(Bl) 



where 



M 
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g + cr s - a. 



P ^-0 2 



g + a a M 16 = a c q 



Mi 3 = a'g 2 



2 a 2 -</> 2 



M 55 = 2a, 



M 33 - g 4 + ( a s - a v ^ - ) q 2 



The vector v q multiplying the matrix M g consists of the real and imaginary parts of the 
modes h q , u q and Sm^; 
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r) P , (««),, (/i g ) r , (&<,)., (5m+) , {6m+) , [5m q ) r , {Sm q ) t 



(B2) 



Appendix C: The minimised values E m j n , A m j n and s 



qrmn 



The minimised values of E, A and s 9 are given by; 



E 



Op0w' + 2 2 tan (0 O ) - ojj/i' 
a p (f)h' — 2 (j p tan (6*0) — <r 2 u' 



(CI) 



or 



a p (pu' — 2 <T P tan (9q) — a 2 h 



2 V 
p 

I a l ~ ^ \ (T p <ph' + 2<j) 2 tan (0 O ) - a 2 u' 



(C2) 
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and 

Sqmin = ^"7j _ Q2Y (• ' 
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